#if !defined(MTX_TYPES_H)
#define MTX_TYPES_H

#include <cassert>
namespace std { struct forward_iterator_tag; }

namespace mtx {

typedef unsigned int	size_type;	
typedef signed int		diff_type;

struct index {
	size_type begin, end;
	diff_type inc;
	explicit index(size_type b = 0, size_type e = 0, diff_type i = 1)
		: begin(b), inc(i), end(e) { }
	size_type size() const
		{ return (end - begin) / inc; }
};

template <class T> class const_slice {
public:
	class const_iterator;
private:
	const T * _t;
	index _i;
public:
	const_slice()
		: _t(0), _i() { }
	const_slice & init(const T * t, index i)
		{ _t = t; _i = i; return * this; }
	const_iterator begin() const
		{ assert(_t); return const_iterator(_t->begin() + _i.begin, _i.inc); }
	const_iterator end() const
		{ assert(_t); return const_iterator(_t->begin() + _i.end, _i.inc); }
	size_type size() const
		{ return _i.size(); }
public:
	class const_iterator {
	private:
		typename T::const_iterator _val;
		diff_type _inc;
	public:
		const_iterator() { }
		const_iterator(const typename T::const_iterator & val, diff_type inc)
			: _val(val), _inc(inc) { }
		bool operator == (const const_iterator & it) const
			{ return _val == it._val; }
		bool operator != (const const_iterator & it) const
			{ return _val != it._val; }
		const_iterator & operator ++ ()
			{ _val += _inc; return * this; }
		const_iterator & operator + (diff_type d)
			{ _val += _inc * d; return * this; }
		const_iterator & operator - (diff_type d)
			{ _val -= _inc * d; return * this; }
		const_iterator & operator -- ()
			{ _val -= _inc; return * this; }
		const double & operator * ()
			{ return * _val;  }
	};
};

template <class T> class slice {
public:
	class iterator;
	class const_iterator;
private:
	T *	_t;
	index _i;
public:
	slice()
		: _t(0), _i() { }
	slice & init(T * t, index i)
		{ _t = t; _i = i; return * this; }
	const_iterator begin() const
		{ assert(_t); return const_iterator(_t->begin() + _i.begin, _i.inc); }
	iterator begin()
		{ assert(_t); return iterator(_t->begin() + _i.begin, _i.inc); }
	const_iterator end() const
		{ assert(_t); return const_iterator(_t->begin() + _i.end, _i.inc); }
	iterator end()
		{ assert(_t); return iterator(_t->begin() + _i.end, _i.inc); }
	size_type size() const
		{ return _i.size(); }
public:
	class iterator {
	public:
		typedef std::forward_iterator_tag	iterator_category;
		typedef double value_type;
		typedef double * pointer;
		typedef double & reference;
		typedef diff_type difference_type;
	private:
		typename T::iterator _val;
		diff_type	_inc;
	public:
		iterator() { }
		iterator(const typename T::iterator & val, diff_type inc)
			: _val(val), _inc(inc) { }
		bool operator == (const iterator & it) const
			{ return _val == it._val; }
		bool operator != (const iterator & it) const
			{ return _val != it._val; }
		iterator & operator ++ ()
			{ _val += _inc; return * this; }
		iterator & operator += (diff_type d)
			{ _val += _inc * d; return * this; }
		iterator operator + (diff_type d)
			{ return iterator(_val + _inc * d, _inc); }
		iterator & operator -- ()
			{ _val -= _inc; return * this; }
		iterator & operator -= (diff_type d)
			{ _val -= _inc * d; return * this; }
		iterator operator - (diff_type d)
			{ return iterator(_val - _inc * d, _inc); }
		double & operator * ()
			{ return * _val;  }
	};
	class const_iterator {
	public:
		typedef std::forward_iterator_tag	iterator_category;
		typedef const double value_type;
		typedef const double * pointer;
		typedef const double & reference;
		typedef diff_type difference_type;
	private:
		typename T::const_iterator _val;
		diff_type _inc;
	public:
		const_iterator() { }
		const_iterator(const typename T::const_iterator & val, diff_type inc)
			: _val(val), _inc(inc) { }
		bool operator == (const const_iterator & it) const
			{ return _val == it._val; }
		bool operator != (const const_iterator & it) const
			{ return _val != it._val; }
		const_iterator & operator ++ ()
			{ _val += _inc; return * this; }
		const_iterator & operator += (diff_type d)
			{ _val += _inc * d; return * this; }
		const_iterator operator + (diff_type d)
			{ return const_iterator(_val + _inc * d, _inc); }
		const_iterator & operator -- ()
			{ _val -= _inc; return * this; }
		const_iterator & operator -= (diff_type d)
			{ _val -= _inc * d; return * this; }
		const_iterator operator - (diff_type d)
			{ return const_iterator(_val - _inc * d, _inc); }
		const double & operator * ()
			{ return * _val;  }
	};
};

template <class L> inline void vassign(L first, L last, double x)
	{ for (L y = first; y != last; ++y) (* y) = x; }

template <class L, class R> inline void assign(L first, L last, R x)
	{ for (L y = first; y != last; ++y, ++x) (* y) = (* x); }

template <class L, class R> inline void rev_assign(L first, L last, R x)
	{ for (L y = first; y != last; ++y, ++x) (* y) = -(* x); }

template <class L, class R> inline void add_assign(L first, L last, R x)
	{ for (L y = first; y != last; ++y, ++x) (* y) += (* x); }

template <class L, class R> inline void sub_assign(L first, L last, R x)
	{ for (L y = first; y != last; ++y, ++x) (* y) -= (* x); }

template <class L, class R> inline void mul_assign(L first, L last, R x)
	{ for (L y = first; y != last; ++y, ++x) (* y) *= (* x); }

template <class L> inline void mul_vassign(L first, L last, double x)
	{ for (L y = first; y != last; ++y) (* y) *= x; }

template <class L> inline void div_vassign(L first, L last, double x)
	{ for (L y = first; y != last; ++y) (* y) /= x; }


template <class T> class vector : public T {
public:
	typedef vector< slice<T> > subvector;
	typedef vector< const_slice<T> > const_subvector;
public:
	vector()
		: T() { }
	vector(double x) : T()
		{ * this = x; }
	inline vector & operator = (double x)
		{ vassign(T::begin(), T::end(), x); return * this; }
	template <class U> inline vector & operator = (const vector<U> & rhs)
	{
		std::copy(rhs.begin(), rhs.end(), T::begin());
		/*assign(T::begin(), T::end(), rhs.begin());*/ return * this;
	}
	template <class U> inline vector & operator += (const vector<U> & rhs)
		{ add_assign(T::begin(), T::end(), rhs.begin()); return * this; }
	template <class U> inline vector & operator -= (const vector<U> & rhs)
		{ sub_assign(T::begin(), T::end(), rhs.begin()); return * this; }
	inline vector & operator *= (double x)
		{ mul_vassign(T::begin(), T::end(), x); return * this; }
	inline vector & operator /= (double x)
		{ div_vassign(T::begin(), T::end(), x); return * this; }
	const double & operator () (size_type i) const
		{ return T::operator () (i); }
	double & operator () (size_type i)
		{ return T::operator () (i); }
	const_subvector operator () (index ix) const
		{ const_subvector v; v.init(this, ix); return v; }
	subvector operator () (index ix)
		{ subvector v; v.init(this, ix); return v; }

};

template <class T>
vector<T> operator + (const vector<T> & x, const vector<T> & y);

template <class T>
vector<T> operator - (const vector<T> & x, const vector<T> & y);

template <class T>
vector<T> operator - (const vector<T> & x);

template <class T>
vector<T> operator * (const vector<T> & v, double x);

template <class T>
vector<T> operator / (const vector<T> & v, double x);

template <class T>
vector<T> operator - (const vector<T> & x)
{
	vector<T> t;
	rev_assign(t.begin(), t.end(), x.begin());
	return t;
}

template <class X, class Y>
double operator * (const vector<X> & x, const vector<Y> & y) {
	assert(x.size() == y.size());
	double buf(0);
	typename X::const_iterator xit = x.begin();
	typename Y::const_iterator yit = y.begin();
	for (; xit != x.end(); ++xit, ++yit)
		buf += (* xit) * (* yit);
	return buf;
}

template <class T>
double abs(const vector<T> & v) {
	return sqrt(v * v);
}

template <class T> class matrix : public T {
public:
	typedef vector< slice<T> > subvector;
	typedef vector< const_slice<T> > const_subvector;
public:
	matrix()
		: T() { }
	matrix(double x) : T()
		{ * this = x; }
	inline matrix & operator = (double x)
		{ vassign(T::begin(), T::end(), x); return * this; }
	template <class U> inline matrix & operator = (const matrix<U> & rhs)
		{ assign(T::begin(), T::end(), rhs.begin()); return * this; }
	template <class U> inline matrix & operator += (const matrix<U> & rhs)
		{ add_assign(T::begin(), T::end(), rhs.begin()); return * this; }
	template <class U> inline matrix & operator -= (const matrix<U> & rhs)
		{ sub_assign(T::begin(), T::end(), rhs.begin()); return * this; }
	inline matrix & operator *= (double x)
		{ mul_vassign(T::begin(), T::end(), x); return * this; }
	inline matrix & operator /= (double x)
		{ div_vassign(T::begin(), T::end(), x); return * this; }
	const double & operator () (size_type i, size_type j) const
		{ return T::operator () (i, j); }
	double & operator () (size_type i, size_type j)
		{ return T::operator () (i, j); }
	const_subvector operator () (size_type i, index ix) const
		{ return subrow(i, ix); }
	subvector operator () (size_type i, index ix)
		{ return subrow(i, ix); }
	const_subvector operator () (index ix, size_type j) const
		{ return subcol(j, ix); }
	subvector operator () (index ix, size_type j)
		{ return subcol(j, ix); }
	const_subvector subrow(size_type i, index ix) const
		{ const_subvector v; v.init(this, index(i + T::rows() * ix.begin, i + T::rows() * ix.end, T::rows() * ix.inc)); return v; }
	subvector subrow(size_type i, index ix)
		{ subvector v; v.init(this, index(i + T::rows() * ix.begin, i + T::rows() * ix.end, T::rows() * ix.inc)); return v; }
	const_subvector subcol(size_type j, index ix) const
		{ const_subvector v; v.init(this, index(j * T::rows() + ix.begin, j * T::rows() + ix.end, ix.inc)); return v; }
	subvector subcol(size_type j, index ix)
		{ subvector v; v.init(this, index(j * T::rows() + ix.begin, j * T::rows() + ix.end, ix.inc)); return v; }
	const_subvector row(size_type i) const
		{ return subrow(i, index(0, T::cols())); }
	subvector row(size_type i)
		{ return subrow(i, index(0, T::cols())); }
	const_subvector col(size_type j) const
		{ return subcol(j, index(0, T::rows())); }
	subvector col(size_type j)
		{ return subcol(j, index(0, T::rows())); }
};			

template <class T>
matrix<T> operator - (const matrix<T> & x)
{
	matrix<T> t;
	rev_assign(t.begin(), t.end(), x.begin());
	return t;
}

template <class L, class R>	void copy(L & lhs, const R & rhs)
	{ assign(lhs.begin(), lhs.end(), rhs.begin()); }

} // mtx

#endif   // MTX_TYPES_H
